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Hanbury-Brown-Twiss interferometry is a technique which yields effective widths (i.e., “HBT 
radii”) of homogeneity regions in the fireballs produced in heavy ion collisions. Because the initial 
conditions of these collisions are stochastically fluctuating, the measured HBT radii also exhibit 
variation on an event-by-event basis. However, HBT measurements have, to date, been performed 
only on an ensemble-averaged basis, due to inherent limitations of finite particle statistics. In this 
paper, we show that experimental measurements to date are best characterized theoretically as 
weighted averages of the event-by-event HBT radii, and we propose a new method for extracting 
experimentally both the arithmetic mean and the variance of the event-by-event distribution of HBT 
radii. We demonstrate the extraction of the mean and variance of this distribution for a particular 
ensemble of numerically generated events, and offer some ideas to extend and generalize the method 
to enable measurement of higher moments of the HBT distribution as well. 

PACS numbers: 25.75.-q, 12.38.Mh, 25.75.Ld, 24.10.Nz 


I. INTRODUCTION 

HBT interferometry relies on two-particle momentum 
correlations to extract information about the spatiotem- 
poral structure of the emitting source in heavy-ion colli¬ 
sions. The technique depends on the detection of pairs 
of identical particles (e.g., pions or kaons), whose quan¬ 
tum statistical correlations convey important informa¬ 
tion about the mean relative separation between the 
points at which the particles were emitted during the 
freeze-out process. Ideally, one would be able to do this 
on an event-by-event basis: if any given collision yielded 
a sufficiently large number of the desired particles in the 
final state, this would allow a measurement of the HBT 
radii, the effective widths of the homogeneity regions in 
the fireball, event by event. Unfortunately, after the to¬ 
tal particle multiplicity (on the order of a few 1000 per 
event) is binned according to particle species, pt, and 
emission angle, not enough pairs remain for a statisti¬ 
cally meaningful, fully three dimensional analysis of the 
correlation function. 

Consequently, experimentalists typically combine large 
numbers (> 10®) of events in order to boost the pair 
statistics, thereby increasing the precision of the result¬ 
ing HBT measurements. The collection of events is re¬ 
ferred to as the ensemble, and the two-particle correlation 
function (from which the HBT radii are experimentally 
extracted) thus contains a non-trivial combination of the 
correlation functions of all of the events in the ensemble. 
An apples-to-apples comparison with theoretical models 
therefore requires, at least in principle, a corresponding 
ensemble averaging on the theoretical side. 

The process of ensemble averaging has historically been 
accounted for at the level of the initial state of the hre- 
ball. In its crudest form, the ensemble of fluctuating 
events is replaced by a single averaged event whose final 
state is computed by hydrodynamically evolving a single 
averaged initial profile. With the recent availability of 
resources to evolve large numbers of collisions with fluc¬ 


tuating initial conditions event by event, the ensemble av¬ 
eraging procedure has been shifted to the emission func¬ 
tion. By performing an ensemble average directly over 
emission functions constructed from the freeze-out sur¬ 
faces of each event in the ensemble, the two particle cor¬ 
relation function can be related to the Fourier transform 
of the ensemble-averaged emission function. However, 
since the experimental correlation function is constructed 
after the final state particles have been emitted from the 
freeze-out surfaces of each respective event, it is more ac¬ 
curate to perform the ensemble-averaging procedure at 
the level of the correlation function itself. This induces 
corrections to the HBT radii extracted from the corre¬ 
lation function which is constructed from the ensemble- 
averaged emission function, and these corrections are sen¬ 
sitive to event-by-event fluctuations encoded in the struc¬ 
ture of the freeze-out surface. Only this last procedure 
invokes a distribution of correlation functions and thus 
a distribution of HBT radii which can be characterized 
by a mean, a variance, and possible higher moments. In 
this paper we will analyze which moment of this distribu¬ 
tion is represented by the experimentally measured HBT 
radii, and what additional measurements could be made 
to access other moments of the HBT radii distribution. 

Numerical studies such as [1] have shown that the mean 
HBT radii extracted from a fluctuating set of correla¬ 
tion functions are almost indistinguishable from the radii 
characterizing the single correlation function obtained 
from an ensemble-averaged emission function. On the 
other hand, recent event-by-event simulations [2] indi¬ 
cate that the HBT radii extracted from individual events^ 
may fluctuate with a typical range of 10-15% (in the 
squared HBT radii) for central collisions and that, if 


^ This is possible in theory since the correlation function does not 
need to be sampled with a finite number of particles but can be 
calculated with infinite precision. 
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these fluctuations are both present in actual heavy-ion 
collisions and experimentally accessible, their scale could 
provide valuable sensitivity, e.g., to different functional 
forms of the T-dependence in vy/s. In this paper, we 
propose a method for extracting the scale of these fluc¬ 
tuations experimentally. A more thorough exploration of 
our method’s theoretical implications is deferred to an¬ 
other work. 

The layout of the remainder of this paper is as fol¬ 
lows. In section II, we introduce the formalism under¬ 
lying both event-by-event and ensemble-averaged HBT 
analyses which retain the correlation function’s intrin¬ 
sic dependence on the pair emission angle, and distin¬ 
guish several different, commonly used methods for the¬ 
oretically computing the ensemble-averaged HBT radii 
which are measured experimentally. We argue that each 
of these methods traces qualitatively, but not with quan¬ 
titative precision, the arithmetic mean of the event-by- 
event distribution of the HBT radii. In section HI, we 
discuss analogous results for HBT analyses in which the 
dependence on the pair emission angle is averaged over, 
and show how this simplification affects the differences 
between the various methods of ensemble averaging. The 
best theoretical representation of the experimentally em¬ 
ployed ensemble averaging process identifies the mea¬ 
sured HBT radii as weighted averages of the event-by- 
event radii. We discuss these weights. With this in 
mind, we show in section V how to estimate the first mo¬ 
ment (i.e., the mean) of an event-by-event distribution in 
terms of linear combinations of such weighted averages. 
Although we focus in this paper on the HBT radii, the 
methods we introduce are quite general and applicable to 
any event-by-event observables. In sections VI and VH, 
we show how to access higher moments of the event-by- 
event distribution by performing repeated averages over 
sub-ensembles and measuring their fluctuations; in sec¬ 
tion VI, we concentrate on estimating the variance, while 
section VH addresses higher moments. Finally, in sec¬ 
tion VHI, we present a proof-of-principle demonstration 
of our method and discuss some subtleties relevant for 
its experimental implementation. A detailed derivation 
of our method (for estimating the variance) is provided 
in Appendix A. 


A. Azimuthally sensitive interferometry for a 
single fluctuating event 


HBT interferometry is founded on the concept of the 
two-particle correlation function, defined for a single 
event by 


C{pi,P2) 


^Pl^P2 d3pid3p2 




d^N A 

^P2 d^p2 I 


( 1 ) 


Here, pi and p 2 represent the 3-momenta of identical par¬ 
ticles (e.g., pions) which have been emitted from the fire¬ 
ball. The correlation function (1) may be interpreted as 
the probability of simultaneously measuring two parti¬ 
cles with momenta pi and p 2 in a single event, divided 
by the probability of measuring the same two particles 
(with the same momenta) independently in two separate 
but identical events. Correlations among the particles 
in the emitted pair manifest themselves as deviations of 
C{pi,p 2 ) from unity. The connection of C{pi,p 2 ) with 
the size of the effective emission region (’’homogeneity 
region”) from which the pairs are emitted is provided 
by the following connection [5] with the single-particle 
Wigner density (or ’’emission function”) of the fireball, 

six,Ky.‘^ 


C(q,K) 


1± 


/ d^xS{x,K)e^<i-^ 
Jd‘^xS{x,K) 


( 2 ) 


where we have introduced the notation = p^ — P 2 , 
AT^ = (pj -|-p2)/2. Eq. (2) holds in the absence of final 
state interactions between the emitted particles and for 
“chaotic” sources that emit the two particles indepen¬ 
dently from each other. The approximation indicated 
by the « sign refers to the replacement of pi,P 2 by K 
in the denominator (the so-called “smoothness approx¬ 
imation” [5]). If the pairs of identical particles used in 
the construction of the numerator and denominator of 
Eq. (I) are bosons (as we consider in this paper), the cor¬ 
relation function itself experiences an enhancement near 
q = 0; this enhancement is usually described by a func¬ 
tional form which is Gaussian in the components of the 
relative momentum q: 


II. AZIMUTHALLY SENSITIVE HBT 
INTERFEROMETRY FOR FLUCTUATING 
SOURCES 

In this section we discuss HBT interferometry that is 
fully differential in the pair momentum AT, in particular 
its azimuthal angle around the beam direction. This 
is known as “azimuthally sensitive HBT interferometry” 
[3, 4]. In Sec. HI we will modify the treatment for az¬ 
imuthally averaged (i.e., <i>if-integrated) measurements. 


Cf.tiq, K) = l + X{K) exp I - ^ Rl^{K)q,qj . 

\ i,j=o,s,l J 

_ ( 3 ) 

This Gaussian parametrization is exact for emission func¬ 
tions with a Gaussian spatial structure and is usually ad¬ 
equate for non-Gaussian sources whose deviations from 
Gaussian structure are generated by additional length 


^ The sign corresponds to using bosons to construct the corre¬ 
lator, whereas the corresponds to using fermions. 
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scales characterizing the source that are very different 
from the source radii. Here, \{K) (the “intercept param¬ 
eter”) encodes information about long-lived resonances 
which decay well outside the reaction zone of the fireball, 
and is used to account for the resulting empirical reduc¬ 
tion in the peak value of C at tf = 0 [5, 6]. We neglect 
the contributions from resonances and set \{K) = 1. The 
sum in the exponent ranges over the coordinates of the 
widely used os/-system, where I (the “longitudinal” direc¬ 
tion) coincides with the beam direction, o (the “outward” 
direction) points in the same direction as Kt, the aver¬ 
age pair momentum projected onto the transverse plane, 
and s (the “sideward” direction) points perpendicular to 
both of these. In terms of these coordinates, Eq. (3) de¬ 
fines the HBT radius parameters whose diagonal 

components may be interpreted as the squares of the ef¬ 
fective widths of the emission regions within the fireball 
responsible for producing particle pairs with average mo¬ 
mentum K. 

One prescription for computing the HBT radii from 
theoretical (hydrodynamical) models on an event-by- 
event basis relies on the Cooper-Frye formula [7] to de¬ 
fine the eventwise emission function, and then uses (2) to 
define the corresponding correlation function. This cor¬ 
relation function can then be fit using (3) to obtain the 
Rfj for that event as fit parameters. Explicitly, in the 
Cooper-Frye algorithm the emission function is defined 
as follows: 

S{x,p) = J^p-d^a{y)S‘^{x-y)fiy,p), (4) 

f{x,p) = fo{x,p) + 5f{x,p) 

_ 1 ^ t , t \ 

+ 2T^{e+VY'^^ ^ 

Here, 6f is the first-order viscous correction to the local 
equilibrium distribution function /o [8, 9], and we assume 
a quadratic dependence on p. tt^v{x) is the viscous pres¬ 
sure tensor, u^{x) is the flow velocity profile along the 
freeze-out surface, and p, T^ec, e, and V are the chem¬ 
ical potential, decoupling temperature, energy density, 
and pressure, respectively, which, for the hydrodynam¬ 
ical simulations used in this work, will all be taken as 
constant along the freeze-out surface by construction. S 
represents the freeze-out hypersurface over which the in¬ 
tegration is performed, and d^a^x) is the outward point¬ 
ing normal vector at the point x on this surface. 

This prescription for generating from theoretical 
models for comparison with experimental data is compu¬ 
tationally intensive: to construct the correlation function 
for a given event according to (2) requires multidimen¬ 
sional integrations over the freeze-out surface of the event 
in question (where S{x,K) itself, in general, requires a 
similar integration if computed according the Cooper- 
Frye prescription (4) and (5)); these integrations, more¬ 
over, must be performed for a sufficiently large number 
of points in q and K that the fit to (3) can be carried 
out with acceptable accuracy and still yield useful results 


for comparison with experiment. To do all of this for a 
large ensemble of events consequently places stringent 
demands on available computational resources. 

Much of this numerical expense can be avoided by 
adopting the following often-used approximation. By as¬ 
suming that the emission function (and, consequently, 
the corresponding correlation function) for each event 
can be described exactly as a Gaussian in x (corresp., 
q), the fit of (2) to (3) becomes an identity, and the 
may be read off directly in terms of integrals over the 
emission functions of the ensemble of events [5, 6]: 


R%{K) = {{ii - - I3ji))s > 


where - (x^)^, fd = K/Ek, and 


{f{x))s = 


f d^x /(x) S{x, K) 
fd^xS(x,K) 


( 6 ) 

(7) 


Although still computationally intensive, the project of 
performing event-by-event HBT analyses has been cast 
into a much more tractable form by the use of (6) and 
(7), which drastically reduce the number of required cal¬ 
culations for computing the from theoretical models. 
Although for realistic (i.e., hydrodynamic) model sources 
there exist well-documented discrepancies between the 
“Gaussian fit” method and “source variances” method, 
these discrepancies will not affect the qualitative results 
that we discuss in this paper, which will be based on 
Eqs. (6) and (7). 


B. Azimuthally sensitive interferometry for 
ensembles of fluctuating events 


While theoretically well-defined on an event-by-event 
basis, single-event HBT interferometry is in general not 
practically possible, as explained in the introduction. For 
this reason, experimentalists typically modify the defini¬ 
tion of the correlation function to include in the numer¬ 
ator and denominator of (1) pairs of pions from multiple 
events; the collection of all events combined in this way is 
known as the ensemble, and the corresponding definition 
of the correlation function is 


Gavg(Pl)P2) — 


E, 


<fN \ 


Pi-^P2 d^pid^p2 j^ 
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Ppi 


E, 


PN \ 
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( 8 ) 


where the (.. Yv notation is shorthand for 


JVe. 




(9) 


i=l 


i.e., an arithmetic average of the quantity X over all 
events in the ensemble. In the construction of the cor¬ 
relator according to Eq. (8), it is important to align the 
events according to some direction defined by the indi¬ 
vidual event, for example, the nth order flow angle 'I'n of 
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charged hadrons. This will result in an azimuthal depen¬ 
dence of the HBT radii relative to that flow angle 
In Appendix B, we show that the radii extracted from 
a Gaussian fit of (8), which we label by can be 

related to the event-by-event HBT radii by 


RhiK) 


(N ^ (i^) y 

_\_ / ev 


( 10 ) 


where N{K) = EK{d^N/(PK) is the Lorentz-invariant 
yield of the particles of interest (in our case, pions) with 
momentum K. Theoretically, there are many different 
ways to generalize (1) to ensembles containing multiple 
events. In the rest of this subsection, we discuss several 
of these alternatives. 


1. Ensemble-averaged initial conditions - single-shot 
hydrodynamics 

One method for extracting from hydrodynamic codes 
a correlation function containing averaged information 
from multiple events involves ensemble-averaging (possi¬ 
bly after proper alignment if the events are deformed) 
the initial entropy (or energy) densities of all events in 
the transverse plane, and using the resulting averaged 
density profile as a set of initial conditions for hydro¬ 
dynamics. At the end of the hydrodynamic evolution, 
Eq. (3) relates the emission function constructed on the 
freeze-out surface for this averaged density prohle to a 
correlation function (by (2)) which is effectively insensi¬ 
tive to the existence of event-by-event fluctuations. We 
refer to this method as ” single-shot hydrodynamics”, and 
we refer to the emission function (resp., correlation func¬ 
tion) so constructed as Rgsh (resp., Cgsh), and denote the 
HBT radii extracted with this method by R?^.. Explicitly, 
we may write (using the shortcut discussed above) 

= {{xi - Pii){xj - , ( 11 ) 

where (.. is defined with Sssh as the weight function. 


2. Ensemble-averaged emission function 


Another common way of computing an ensemble- 
averaged correlation function consists of averaging the 
emission functions after the event-by-event hydrody¬ 
namic evolution of many events with fluctuating ini¬ 
tial conditions. In particular, we define S{x,K) = 
{S(x, K))^^ and, in analogy with (2), introduce 

Jd^xSix_,K)E<^- " 

f d'^x S{x, K) 


C{q,K) = l + 


corresponding “single-shot hydrodynamics” definitions 
that we discussed above, in that the quantum fluctua¬ 
tions in the initial state are allowed to modify the hydro- 
dynamic evolution event by event before being averaged 
over at the end. Since the hydrodynamic evolution is 
nonlinear, S{x,K) ^ Ssshix, K). In this work, we will 
refer to this method as the ’’average emission function” 
method, and dehne 

Ry (R) = ((ii-At)(ij-/?jt))s . (13) 


3. Ensemble-averaged correlation function 

Of the available theoretical techniques for treating the 
ensemble-averaging process, the “single-shot hydrody¬ 
namics” and “average emission function” methods have 
historically enjoyed the greatest popularity. However, as 
Eq. (8) shows, the way to correctly reproduce the exper¬ 
imental process of performing the ensemble average is by 
first Fourier transforming the emission function and then 
averaging over events: 

C'avg(pi,P2) = {C)^^{q,K) 

(l/ d‘^xS{x,K)P’^-^f) 

= 1+^- ^(14) 

1/ d^x 5'(x, K)\ 

1 / d‘^x5S{x,K)E‘‘t'^ 

\j'd‘^xS{x,K)\^ 

where 5S{x, K) = S{x, K) — S{x, K). Using the shortcut 
(6) we can then write the corresponding HBT radii R^^-) 

extracted from a Gaussian fit of (C')g^ (QjK) as follows 
(see Eq. (10)): 

^(b>(^) = , (15) 

where we suppressed the AT-dependence on the right- 
hand side. We refer to this way of computing the HBT 
radii as the “average correlation function” method. As 
an additional check on this result, we note that if we 
neglect event-by-event fluctuations entirely by setting 
SS{x,K) = 0, then S{x,K) = S{x,K), the final term in 
(14) vanishes, and (15) reduces to (13), as expected. The 
correct theoretical definition of the ensemble-averaged 
HBT radii may therefore be thought of as a simple 
(weighted) average of the HBT radii for each fluctuating 
event, each scaled by a factor which accounts explicitly 
for final-state multiplicity fluctuations, bin by bin in K, 
from event to event. 



= Ciq,K) + 


4- Direct ensemble average 


These definitions of S and C (and the corresponding A more direct route skips the construction of the cor- 
radii, which we denote Ufj) have an advantage over the relation function entirely, and simply averages the radii 
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(6), computed from S{x,K) event-by-event, directly: 


corresponding $/f-averaged single-particle spectra: 


'A)JK) 


1 


11(4)“’(A'). 




(16) 


where denotes the HBT radii of the fcth fluctu¬ 

ating event. (Of course, this could also be done if the 

radii were extracted from a Gaussian fit to the 

correlation function for event k using Eq. (3) with the cor¬ 
responding emission function S^^\x,K). Here, however, 

(yR^j) will be computed via the shortcut (6), with S re¬ 
placed by This “direct ensemble average” is clearly 

equivalent to the arithmetic mean of the event-by-event 
radii, without the additional multiplicity weight in (10). 
Hereafter, we will drop the subscript ’ev’ in the interest of 
notational simplicity, whenever doing so is unambiguous. 


III. AZIMUTHALLY AVERAGED HBT 
INTERFEROMETRY FOR FLUCTUATING 
SOURCES 


A. Azimuthally averaged HBT interferometry for a 
single event 


Azimuthally sensitive HBT interferometry relies fun¬ 
damentally on the construction of the two-particle cor¬ 
relation function given in Eq. (1). Again we first study 
the situation for a single event. Taking (1) as a starting 
point, there are at least two distinct ways of obtaining 
azimuthally averaged HBT radii. The first is to construct 
the full azimuthally dependent correlation function, ob¬ 
taining the azimuthally sensitive HBT radii by fitting 
(1) to the form (3), and then average explicitly over the 
residual dependence on ^k- 

R-ij,oi^r) = d^KRij iKT,^K) 

= {rI{Kt,^k))^^. (17) 


The second way to obtain azimuthally averaged HBT 
radii is to perform the average at the level of the correla¬ 
tion function (1) before fitting HBT radii to it, instead of 
averaging the HBT radii after fitting the correlation func¬ 
tion event by event. Since the correlator is constructed as 
the ratio of two experimental quantities which are mea¬ 
sured on a bin-by-bin basis, binning only on Kt without 
binning on avoids reducing the number of available 
particle pairs per bin by a factor of the number of bins in 
^K- Theoretically, a correlator constructed in this way 
should be written as the ratio of the $/f-averaged two- 
particle cross section divided by an uncorrelated back¬ 
ground which is constructed by taking a product of the 





( 18 ) 

In terms of the emission function S{x, K), this correlator 
may be written as 


(C)^^ {q,K) 


(\J d‘^xS{x,K)\^^^ 
{Jd^xS{x,K))^^ 


(19) 


x\l + 


<J|/ d‘^xe^‘^-^Six,K)f'j^ 
(^\J d'^xSix,K)\^'j^ 


As before, we consider fitting this expression to the form 
(3) and extracting the Rfj as fit parameters which de¬ 
pend only on Kt- However, since the factor inside the 
parentheses in (19) tends to 2 in the limit that g —)■ 0, we 
must include an overall factor when fitting the correlator: 


(C')$^(9,l?) = C'o |l+exp| - q^qjRljiKr) 

\ \ i,j=o,s,l 

where 


( 20 ) 


{NHKT,<i>K))^ 

{N{KT,^K))'i^ 


( 21 ) 


with N{Kt,^k) = f d'^xS{x, K) as before. The az¬ 
imuthally averaged R^j^Kt) are proportional to the cur¬ 
vature of the correlator at the origin, and one can 
show that this leads (with the functional dependence of 
S{x, K) suppressed) to 


RI{Kt) 


{J d^x ixi-l3^t){xj-l3jt)S)^^ 

if d^x {xi-p,t)S)^^ (f d^x {xj-Pjt)S)^ 

(fd^xS)l^ 

(N'^(Kt, ^K)Rfj(KT, ^k))^ 




( 22 ) 


Eor Gaussian sources with azimuthally symmetric par¬ 
ticle emission, the definition (22) is equivalent to (17), 
i.e., Rfj(KT) = RfjQiKx), since the factor N‘^{Kt,^k) 
is -independent and thus drops out from the ratio in 
(22). For events with a significant azimuthal asymme¬ 
try in pair production, however, these two methods of 
azimuthal averaging can yield substantially different re¬ 
sults. 

Note that the prefactor outside the parentheses in (19) 
is independent of q and thus also modulates the correla¬ 
tion function at |(^ —)■ oo. It does not affect the extrac¬ 
tion of the HBT radii if, as is often done in experiment, 
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the correlation function is normalized by hand to 1 at 
1^ oo. 


B. Azimuthally averaged HBT interferometry for 
ensembles of flnctuating events 


1. Azimuthally averaged HBT radii Ft^j Q{KT) 


In the previous subsection, we introduced two different 
ways of defining the azimuthally averaged HBT radii for 
a single event. Analogously, for an ensemble of events, 
there are two different ways to average over the ^k- 
dependence of each of the ensemble averaging methods 
defined in Eqs. (11), (13), (15) and (16). We can either 
first perform an azimuthally sensitive HBT analysis, ex¬ 
tract the <l>/f-independent HBT radii, and average these 
over ^K, or perform the ‘h/^-average already at the level 
of constructing the correlator (see Eq. (18)) and then 
extract <l>/f-independent radii from the azimuthally aver¬ 
aged correlator. The first procedure requires higher event 
statistics, and is therefore experimentally more difficult. 
Still, it is of conceptual interest and will thus be studied 
in this subsection. The second method is experimentally 
preferred and will be discussed in the following subsec¬ 
tion. The “direct ensemble average” of the azimuthally 
averaged HBT radii (17) is given by 

(4,0> ^ (4-,o(^t)L ^ {{Rl (Kt, • 

"(23) 

Since the azimuthal average commutes with the arith¬ 
metic average over events, (Rij) is also the $ic-average 
of the ensemble-averaged azimuthally symmetric ra¬ 
dius (16). Similarly constructed azimuthal averages of 
Eqs. (11), (13) and (15) define their <l>;f-independent 
parts, with 




(24) 

4.o(^) ^ 


(25) 



(26) 


2. HBT radii B^j{KT) from azimuthally averaged 
correlators 


In the same way that our treatment of azimuthally 
sensitive correlation functions in Sec. H admitted several 
different strategies for generalizing to ensemble-averaged 
correlators, the correlator (18) possesses several analo¬ 
gous ensemble-averaged generalizations. 


We define the “direct ensemble average” of the az¬ 
imuthally averaged HBT radii by 

(4> s 

{N'^{Kt, ^k)R%{Kt, $ic))^ 

in accordance with (22). In the cases of the “single-shot 
hydrodynamics” and “averaged emission function” meth¬ 
ods, both approaches can be characterized by a single 
emission function (either Sssh or S), and therefore imply 
the following ensemble averaged generalizations of (22): 

^k)RIj{Kt, ^k) 



RUKt) ^ 




and 


RI{Kt) = 


(^{N{Kt, ^K))ly RljiKr, 
\ / 


(29) 


where NsshiKr, ^k) and {N)^^ {Kt, ^k) are defined in 
an obvious way. 

An equally straightforward, but slightly more tedious 
derivation shows that the “average correlation function” 
method discussed in Sec. HB 3, which simulates the pro¬ 
cedure applied in experimental analyses [10-12], can be 
adapted to the azimuthally averaged case by defining 




((iV2(AT,<i>K))ev)^^ 

(rHKt) 


Note that event-by-event fluctuations of the spectrum 
N{Kt,^k) only enter in this last method which repro¬ 
duces the experimental procedure. 


IV. COMPARISON OF THEORETICAL 
ENSEMBLE AVERAGING PROGEDURES 

In general, each of the ensemble averaging methods dis¬ 
cussed above yields different results. While experimen¬ 
tally only the “averaged correlation function” method 
is available, leading to the two possible ways (24) and 
(30) to measure the HBT radii, the “single-shot hydro¬ 
dynamics” (Eqs. (25) and (28)) and “ensemble-averaged 
emission function” methods (Eqs. (26) and (29)) can be 
used in theoretical studies and offer significant numer¬ 
ical advantages. It is therefore of interest to evaluate 
the significance of the differences between the different 
ensemble-averaging prescriptions. We begin by compar¬ 
ing in Fig. 1 the HBT radii {Rjj^), Rfij)^^ 
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R^. p (defined in Sec. IIIB 1, Eqs. (23) - (26)) for each of 
the four prescriptions, applied to a typical hydrodynamic 
analysis using the iEBE-VISHNU package [13]. Here, we 
consider 200 A GeV Au+Au collisions at 0-10% central¬ 
ity, using the MC-Glauber model with p-l-p multiplicity 
fluctuations to compute the fluctuating initial entropy 
density profiles in the transverse plane, evolving them 
with boost-invariant hydrodynamics (with 77/5 = 0, i.e., 
assuming ideal fluid behavior) to simulate the evolution 
of the fireball.^ We use Ng^ = 5000 which is large enough 
such that the observed variance of the HBT radii is dom¬ 
inated by event-by-event fluctuations, and fluctuations 
from finite sampling statistics can be neglected. We ter¬ 
minate the hydrodynamical evolution along a freeze-out 
surface of constant temperature Tdec = 120 MeV and use 
the Cooper-Frye algorithm to compute the charged parti¬ 
cle yields. Additional details of our analysis, as well as a 
more systematic discussion of the effects of shear viscos¬ 
ity on HBT analyses, are described in [2]. Fig. 1 shows 
that single-shot hydrodynamics, Rf^ p, (red dash-dotted 
line) leads to the least reliable theoretical estimates for 
the directly ensemble-averaged HBT radii This 

reflects the strongly non-linear hydrodynamic response 
to event-by-event fluctuations in the initial density profile 
which single-shot hydrodynamics does not capture. Both 
the “average emission function” R^j^Q (blue dashed line) 
and “average correlation function” p (green dotted 
line) methods yield results that are in much better agree¬ 
ment with the direct ensemble average (.Rfj,o)ev (black 
solid line), although they also tend to deviate from it at 
Kt > 1 GeV. 

In Fig. 2 we show what happens to these four types of 
azimuthally and ensemble-averaged radii if the azimuthal 
average is not performed at the end of the HBT analysis 
on the level of the HBT radii, but instead on the level of 
constructing the correlation function, before extracting 
the HBT radii, as discussed in Sec. IIIB2, Eqs. (27) - 
(30). We focus our attention on the black solid and green 
dotted curves, representing the algebraic mean and ex¬ 
perimentally determined average radii, {R-fj) and 
respectively. We see that using the appropriate prescrip¬ 
tion for d)if-averaging that applies to each method signif¬ 
icantly improves the agreement between the experimen¬ 
tally accessible HBT radii Rfij'^ and the theoretically in¬ 
teresting algebraic means (Rij), especially at large Kt, 
when compared to the radii defined via Eqs. (27)-(30) 
that were shown in Fig. 1. 

After this discussion of the first moment (i.e., the 


® Taking rj/sf^O tends to suppress the effects of event-by-event 
fluctuations, leading to much smaller discrepancies between the 
various methods of ensemble averaging [2]. The corresponding 
discrepancies amongst the radii derived from ideal hydrodynam¬ 
ics (shown in the figures) therefore represent an upper limit on 
the extent to which the different methods of ensemble averaging 
may disagree with one another for arbitrary values of r]/s. 



FIG. 1: The radius parameters extracted according to the four 
ensemble averaging methods described in Eqs. (23) - (26) (left 
column), as well as their percentage deviation from the direct 
ensemble average (i?fj_o(AT))g^ (right column), as functions 
of Kt- The shaded bands represent the standard uncertainty 
of the direct ensemble average resulting from event-by-event 
fluctuations. 
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FIG. 2: Similar to Fig. 1, but using the definitions (27)-(30) 
for the azimuthally and ensemble-averaged HBT radii. 


mean) of the event-by-event distribution of HBT radii, 
we now proceed to a discussion of the moments of this 
distribution. 









































V. ESTIMATING THE DIRECT ENSEMBLE 
AVERAGE OF THE HBT RADII 


In the preceding sections, we saw that the ensemble of 
events is characterized by an event-by-event distribution 
of HBT radii , and that the experimentally extracted 
ensemble-averaged HBT radii track, but do not exactly 
reproduce the mean value of this distribution. In this 
section, we describe a way of experimentally estimating 
this mean (i.e., the direct ensemble average us¬ 

ing only the experimentally accessible weighted sample 
averages from this distribution. We present this method 
as for a general observable O that can be defined on an 
event-by-event basis, although in this paper we will even¬ 
tually restrict its application to the HBT radii. We will 
refer to this method as “mean estimation.” 

For a general physical observable O that is dehned 
on an event-by-event basis and fluctuates from event to 
event, we distinguish between two types of distributions: 
(i) the underlying physical probability distribution V{0) 
which, for continuous O, is in general a continuous func¬ 
tion, and (ii) the discrete distribution Vn{0) that de¬ 
scribes the distribution of values Oi,... ,On obtained 
in n measurements of the observable O. As n —>■ oo, 
the discrete distribution Vn{0), properly normalized, ap¬ 
proaches the physical distribution 'P{0). We will refer to 
Vn{0) as the “measured distribution of O” or the “en¬ 
semble distribution of O” in our ensemble of n measured 
events, while V{0) will be called the “true” or “physical” 
distribution of O. 

For an arbitrary distribution Q{0), we also dehne for 
later use the associated distribution QniO) as the distri¬ 
bution of sample means of observable O of size n sampled 
from a physical distribution Q{0). Thus, the distribu¬ 
tion of sample means of size n from the true physical 
distribution V{0) above is denoted by Pn{0). For a 
measured distribution Vn{0) of size N we denote the 
analogous distribution of sample means of size n < N 
by pN,n{C>). For sufficiently large ensembles (i.e., as 
N —>• oo), pN,n{0) converges to Pn{0), in the same way 
that Vn{0) converges to V{0). Our goal will be to esti¬ 
mate moments of the ensemble distribution Vn{0) (and 
thereby, the physical distribution V{0)) by studying the 
statistical moments of pN,n{0) in repeated sets of sam¬ 
plings of size n, in the limit of sufficiently large N. 

We now show how to estimate the direct ensemble av¬ 
erage {O) of an observable O. We reiterate that this 
quantity, in the context of HBT interferometry, cannot 
be directly measured experimentally, since single heavy- 
ion collisions yield total multiplicities which are too small 
for extracting meaningful estimates of the radii event by 
event. The fundamental observables available from ex¬ 
perimental HBT analyses are therefore only multiplicity- 
weighted averages of the HBT radii as dehned in Eq. (10). 
We thus consider weighted averages of the observable O 


of the form 

N 

{wO)^ = Y,wi^'^0^, (31) 

k=l 


where k is an index running over all N events, and the 
weights are subject to the following normalization condi¬ 
tion: 


N 

E“-r=i- (32) 

k^l 

For instance, Eq. (10) may be obtained by taking N = 
iVev, O = R^iK), and = Ni{K)/j:^-N^,{K). 

We now show how to construct estimates for the mo¬ 
ments of the event-by-event distribution V{0) in terms 
of expressions of the form (31). 

In order to estimate the direct ensemble average 

1 ^ 

fc=l 


of O in terms of the weighted average {wO )in Eq. (31) 
we must hnd a way to correct for the weights that 
are an unavoidable part of the experimental measure¬ 
ment. Eor the purpose of this paper we will assume that 
the weights are measurable event by event (such 

as the multiplicity weights above). We will also assume 
that for every value w of the weight our measured sample 
contains many events with weights close to w. 

As a first step, let us sort the events by increasing 
weight Next, we create rih bins and hll each bin 

with the same number n = N/nb in order of increasing 
weight \ We label these bins by (£), £ = 1,... ,nt,. 
Eor each bin {£) we construct the weighted average 


- E -’rok 


keii) 


with the modified weights 


(n) _ 

= 


W 


(AT) 


E 


fee(^) ^k 


W 


(34) 


(35) 


is a weighted average of the type that can be 
measured experimentally (such as the HBT radii (10)). 
If the number n of events in each bin is large enough, 
this weighted average will be known with good statistical 
precision, i.e. it will closely approximate the correspond¬ 
ing weighted average of the underlying physical distribu¬ 
tion 'P{0). If n, while being large, is much smaller than 
the total number of events N, the monotonically increas¬ 
ing modified weights will not show much variation^ 


^ If one of the bins happens to contain events from a stretch in the 
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within each bin {£) and will all have approximately the The variance of the ensemble distribution of O is de¬ 
same magnitude 1/n: fined by® 



1 

n 


1U 

~N 


(36) 


We can now arithmetically average the weighted bin av¬ 
erages to find 


ni 

e=i 

kG{i) 






1 ^ 

= — VOfe 

k=l 



^ {0)n 

(37) 


The approximation in the second line above becomes 
exact as N,nb —)■ oo. The ensemble average {O)^ is 
the mean of the ensemble distribution Vn{0) mentioned 
above. For large N it approaches the true mean {O) of 
the observable: 


1 ^ 

al^N ^ Var [Vm{0)] ^ ^ (ol - {0)%) , (39) 

i—k 

where {O) is the ensemble mean of O defined in (37). 
The variance of the physical distribution of O is then 
related to ctq ^ by 

cr^ = lim (40) 

To estimate ctq we assume a very large ensemble of 
N events and consider the process of randomly splitting 
this ensemble into Ub bins of n = N/nt events each, com¬ 
puting the 5u6-ensemble average (0)„ of O for each bin. 
For many different repetitions of this process, the entire 
collection of sub-ensemble averages obtained yields the 
distribution of sample averages of O over sub-ensembles 
of size n from an ensemble containing N events, intro¬ 
duced in Sec. V and denoted by pN,n{0). The central 
limit theorem guarantees that the variance of pN,n{0) is 
proportional to the variance of the ensemble distribution 
of o,rN{oy. 


lim (O)^ = ^ = (O). (38) 

N—¥oo N.rih—^oc) Tlh ‘ ^ 

In Sec. VIII we will show numerical results for a toy ex¬ 
ample where taking nt > 10 in (37) is sufficient to obtain 
an estimate of (O)^ with statistical error of less than 
1 %. 


Ub 


VI. EXTRACTING THE SCALE OF 
FLUCTUATIONS IN THE HBT RADII 


A. Estimating the variance 


We now proceed to show how to estimate the variance 
of the event-by-event distribution of an observable O.® 


(n) 

ordered list where ’ rises abruptly, and thus the condition 
of small variation of the weights inside that bin is violated, this 
bin may be thrown away. In essence, the method described in 
this section corresponds to binning the total ensemble of events 
into bins with approximately the same weight (multiplicity) such 
that within the bin weight (multiplicity) fluctuations can be ne¬ 
glected, and then averaging the bin averages over all bins. If 
total statistics in the ensemble is sufficient, this last average over 
weight bins need not be performed, and this would allow us in 
our case to study whether events with different multiplicities have 
different mean HBT radii. 

^ Note that our use of the term “variance estimation” in this paper 
differs from its more common usage in the field of statistics to 
refer to the general set of techniques for gauging the precision of 
estimators derived from sample data. Although based on simi¬ 
lar principles, the method of variance estimation presented here 
bears only a superficial resemblance to this set of techniques. For 
further discussion see, for instance, [14]. 


Var [PnAO)] cx Var [Vn{0)] . (41) 

In the limit iV —)■ oo, the proportionality factor is simply 
I/n; for finite N, an additional finite population correc¬ 
tion factor must be included (see Sec. VII). 

To formulate an explicit estimate for the variance of the 
ensemble distribution of O, let M be the total number of 
times the ensemble of size N is split into Ub bins of size n, 
and let (O)^, ^ represent the sub-ensemble average of O in 
the £i\i bin of the fcth binning iteration (we suppress the 
dependence of (O)^, ^ on the bin size n to reduce clutter). 
Then we can show (see appendix A) that 

= Mubiub - I) 

is a variance estimator that converges to CTq ^ as M 
approaches the maximal number of different binnings 
M^b,x{N,nb) defined in (A2). According to Eq. (40)), 
CTq ^ itself converges to the variance of the underly¬ 
ing physical distribution V{0) as V —)• oo. Again, we 
suppress the dependence of CTq ^ gg,. on Ub and M for 
clarity. For the particular case of O = (and defining 
afj = 0-^2 ), we have 

AT M rib 

^ Mubiub - 1 ) ^ ^ “ (^Dn) 

° ’ k=i 1=1 

(43) 


® The factor of 1/(A’—1) is used in this expression in place of 1/A^ 
in order to render it an unbiased estimator [15] of the variance 
of the physical distribution of O. 
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Each sub-ensemble average {0)j, ^ may be estimated by 
subdividing the corresponding bin of size n into h\, = njh 
sub-bins of size h and employing the method of mean 
estimation presented in the previous section. Thus, we 
are able to estimate the variance (or, equivalently, the 
second central moment) of the ensemble distribution of 
an observable O by computing the distribution of sub¬ 
ensemble averages (i.e., first moments) from the same dis¬ 
tribution. Furthermore, the equivalence between CTq ^ 
and CTq ^ holds exactly in the limit M —>• M^a^{N,nb). 
However, even for moderate N, Mmax = fV!/(n!)"*’ is 
huge. For example, iV = 10 and nt = 5 yields M^ax ~ 
10®; for N = 100 and ni, = 2, Mmax ~ 10^®. So in prac¬ 
tice the limit M —)■ Mmax is out of reach. Fortunately, 
we will see that ^ gg^ « CTq jy to a very good approxi¬ 
mation, even for M <C Mmax- 

Rather than summing over all possible distinct ways 
of sorting N events into ni, bins, one can thus evaluate 
(42) by summing only over a sufficiently large number 
M of subdivisions such that CTq ^ gg^ converges to a fixed 
value within a given tolerance. In this way, the process 
of evaluating the righthand side of (42) can be performed 
cumulatively. Furthermore, each iteration of the sorting 
and bin-averaging procedures can be performed indepen¬ 
dently of all the others, implying that our method can 
be easily parallelized for numerical computation. 


B. Constructing the covariance matrix 


generalization of (42) to three dimensions.^ 


VII. GENERALIZATION TO HIGHER 
MOMENTS 

We can generalize the combination of the methods in¬ 
troduced in sections V and VI to permit the extraction of 
higher moments of V{0). This method, which we refer 
to as “higher moment estimation,” assumes the ability 
to estimate or calculate the direct ensemble average, and 
therefore relies on the method of mean estimation already 
discussed. We first recast our estimate for the variance 
(43) in terms of the second central moment of Vn{0)'.^ 



= ^2,N,est- (46) 


Here, we have introduced the notation Mk, k > 2, for 
the fcth central moment of the measured O distribution 
Vn{C>), defined by 

1 N 

■ (47) 

^ k=l 


The method discussed in subsection VIA is readily 
extended to incorporate event-by-event fluctuations of 
multiple observables, and the correlations between them. 
These correlations may be quantified by the covariance 
matrix between the observables of interest, and each el¬ 
ement of this matrix may be estimated using a straight¬ 
forward generalization of (42): 


Cov(Oi,02)^_g,,^ 
M rib 

k=l 1 = 1 


N 

nb{nb - l)M 


(Oi)^) ((O2),, 


(44) 



For Oi = and O 2 = Rj’j>, this expression becomes 


Cov 

M rib 

x (<««>».,-<-Rt>„) ((««•>« 



For the most general case in three dimensions, this leads 
to a 6 X 6 covariance matrix for the full set of Rfi. When 
i' = % and = j, this expression reduces to a simple 


The r.h.s. of the (46) consists of two factors: the first is 
a correction factor which accounts for the finite number 
of events N in the ensemble, while the second factor is 
the second central moment of the distribution of bin- 
averages of size n = Njnb, sampled from the ensemble. 
{0)j ^ is a random variable o distributed according to the 
distribution VN,n defined in Sec. V. Defining additionally 
Mi^n = {O )we can thus write 



^ For the discussions and derivations presented so far, we have 
written our results in terms of fluctuating radii Often in the 
literature, however, the reported quantities are not simply 
but Rij, so that one has a choice whether to report properties 
of a distribution of the squared radii or the radii themselves. In 
this paper we will adopt the former approach: the results we 
present are the moments of the event-by-event distribution of 
the fluctuating R^j. 

^ For a smooth distribution, such as the true HBT distribution un¬ 
derlying all HBT measurements, the variance is identically equal 
to the second central moment. For the ensemble of measured 
events, however, the sample size N is large but finite, and the 
variance differs from the second central moment by a factor of 

(V) [151- 
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where (.. .)^ „ denotes the expectation value with respect 
to the distribution VN,n- 

The extension to higher order moments can be found 
in textbooks (e.g., [16]). For the third and fourth order 
moments (related to the skewness and the kurtosis) one 
finds 


\ / N,n 

\ / N,n 


iN-n){N-2n) 


(49) 


n3{N-l){N-2){N-3) 

X [{N^-6nN+N+6n^) M 4 ^Pf 
+3N{n-l){N-n-l) MIj^] , 


which can be solved for and Using our ear¬ 

lier notation we can thus write our estimates for the mo¬ 
ments M 3 TV and M 4 TV as 


Ma^TV.est = 


X 


'n‘^{N-l){N-2) 

^ {N-n){N-2n) 

^ ^ M niy 

EE(<< 3 >w 


UbM 


k=i e=i 



(51) 


M 4 , 


1 


N.est — 


n^N-l){N-2){N-3)\ / 1 


N‘^—6nN+N+6n‘^ \ V N—n 

-3N{n-l){N-n-l) 


M rif, 


UbM 


(^)Af 


niN-l)Y 1 


fc=i e=i 
M rib 


N—n ) \ UhM 


k^i ^=1 




(52) 


In terms of the central moments of a given distribution, 
we can also define its skewness and excess kurtosis: 


0 — ^ 3 ,N 

H 3 ,N — T/9 ’ 

(53) 

0 M4_TV .3 

(54) 


We will find these definitions useful when we demonstrate 
the validity of our method below. Eqs. (53) and (54) 
may also be used to obtain /33,AT,est and /34^TV,est—3, with 
each occurrence of replaced by Mk,N,est as defined 

above. 


VIII. RESULTS AND DISCUSSION 


hydrodynamically evolved fireballs with fluctuating ini¬ 
tial conditions. 

We now illustrate our procedure for estimating the 
arithmetic average from weighted sub-averages. To do 
this, we consider a joint binormal distribution of two ran¬ 
dom variables, which we label X (the observable of in¬ 
terest) and w (the unnormalized weight attached to each 
measurement): 


PYx,iJ'w,<yx,crw,P;X,w) = 
p™ r _ ]_ ( I (x-ux) 

exp 1^ 4_p2 2al + 2<t^ 


(55) 


p{w-fl4n){X-flx) 

<7w<7x 


27 rA/l - p^crx<Jw 


We treat these two variables as governed by a joint dis¬ 
tribution with a non-trivial covariance matrix. 

For specificity, we take the following set of parameters: 


A. Estimating the arithmetic mean 

In this section we present several proof-of-principle 
demonstrations of the various methods discussed in this 
paper. In the present subsection, we illustrate the 
method presented in Sec. V to estimate the arithmetic 
average from particular linear combinations of weighted 
sub-averages. In Sec. VIIIB, we discuss our method for 
estimating higher moments of a distribution by fleshing 
out the sampling distribution of the arithmetic average 
for the same distribution. Finally, in Sec. VIIIC, we 
combine these two methods and use them in a more re¬ 
alistic scenario to estimate the relative variance of event- 
by-event fluctuations in the HBT radii from a sample of 


Hx = 10, ^^ = 7 (56) 

ax = 3,a„ = l/5 (57) 

p = 0.0,0.2,0.4,0.6,0.8,0.999 (58) 

The last variable, p, controls the correlation between the 
stochastic variables X and w. 

We sample from this distribution N = 10,000 

observation-weight pairs {Xi,Wi). Using the procedure 
discussed in Sec. V, we repeatedly distribute these N 
events randomly into rib bins of size h = N/fib and es¬ 
timate the arithmetic mean {X)j^ from Eq. (37). The 
results are plotted in Fig. 3 as a function of the number 
of bins hf,, for several different values of the correlation 
coefficient p. We compute {X)j^ from the lefthand side 
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FIG. 3: The estimated mean (-Y)jvest’ compared with the 
exact mean {X)j^, as a function of hf, and p. We see that de¬ 
creasing the strength of the correlation p leads to a reduction 
in the overall relative difference between (-^)jvest (^)jv> 
since total decorrelation requires that the weighted average 
factorizes and reduces to the arithmetic mean. We note ad¬ 
ditionally that, for hb ^ 10, our method of estimating the 
arithmetic mean is accurate to better than 1%, even with 
nearly total correlation between the Wi and the Xi. 


these fluctuations imposes a lower bound on the error in 
our method of mean estimation which in Fig. 3 is of order 
10“'*. We have found empirically that this lower bound 
on the error scales approximately as 0{\/\/Nhb). The 
accuracy of the method therefore increases as the num¬ 
ber of events in the ensemble N and the number of bins 
fib employed in the method are increased. 


B. Estimating higher moments 

To demonstrate the robustness of our method for ex¬ 
tracting the arithmetic mean and using it to estimate 
higher moments of a distribution, we take a random sam¬ 
ple of iV = 10,000 observations from a skew normal distri¬ 
bution characterized by the following probability density 
function: 


P{fj,,a,a;X) 




TTCr 


2,t2 


^1 -f erf 


/ a{X - fj.) 

V 


(59) 

where /r is a location parameter related to the mean of 
the distribution, a characterizes the width of the distri¬ 
bution, a controls the skewness of the distribution, and 


of Eq. (37) and compare it with the exact mean {X)j^ of 
the N sampled values. 

Figure Fig. 3 demonstrates clearly that for fib ~ 10 we 
can approximate the arithmetic average {X)j^ to better 
than 1%, as noted earlier. All curves eventually drop to 
zero when hb = N = 10,000, since this corresponds to 
placing each event in its own bin, so that the lefthand 
side of (37) reduces to the arithmetic average {X)j^. In 
between, i.e. for 10 < rfh < 10,000, we note that the 
plotted ratio starts to vary irregularly as a function of Ub 
once it reaches a level around 10““*, reflecting the funda¬ 
mental granularity of the measured sample of N “events” 
as explained below. 

In Sec. Ill, we noted that bins which contained 
abruptly changing weights after ordering should be dis¬ 
carded; here, however, we have included all bins for 
the sake of simplicity, regardless of how significantly the 
weights vary within a bin. Fig. 3 demonstrates that this 
choice does not impede our ability to reliably estimate 
the mean {X)j^. 

In the particular case of p = 0 and N ^ oo (keep¬ 
ing hb finite), each weighted bin-average {wOy on the 
lefthand side of (37) factorizes into a product of the av¬ 
erage weight and average observable in each bin. Each 
weighted bin-average thus reduces to the corresponding 
arithmetic bin-average, causing the normalized difference 
between {X) and {X)^, which is plotted in Fig. 3, to 
vanish. The fact that in Fig. 3 the curve for p = 0 does 
not collapse identically to zero can be understood as the 
result of event-by-event fluctuations of the event weights 
about their average values in each bin. The presence of 


erf(z) = —— / dte~* . (60) 

Jo 

For a = 0, iJ, and a correspond to the mean and stan¬ 
dard deviation of P, respectively. In order to illustrate 
the generality of our result, our choice of parameters is 
somewhat arbitrary and entirely without physical moti¬ 
vation: we take p = 17, cr = 3, and a = 10. The proper- 



FIG. 4: The estimated variance compared with the 

exact variance a%, for N = 10,000, as a function of nj, and 

M. 
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ties of this distribution are, for our choice of parameters: 


A^true — 


/ oo 

dx cr, a; x) 

-OO 


= /r + aa 


7r(l + a^) 

2 

/ dx {x — fi) P{^,a,a;x) 

J —OO 


19.3818 (61) 


= ct " 1 - 




TT (1 + a^) 


3.32715 (62) 


/^S.true — 


/^4,true 3 — 


Afs^true 


^true 


■\/2(4 — 7r)a3 


(tt + (tt — 2)a^) 

Mr^true g 

4 "3 

f^true 

8(7r — 3 )q;^ 

(tt + (tt — 2 )a^)^ 


^ « 0.955557 (63) 


0.823244 (64) 


The corresponding statistics for our specific measured 
sample of iV = 10,000 events were 


{X)^ = 19.3588 (65) 

a% = 3.22043 (66) 

= 0.966647 (67) 

^4.Ar- 3 = 1.00149 (68) 


Using the methods defined by (42), (46) and (51)-(54) to 
estimate these sample statistics, for different values of Ub 
and M, we can plot the convergence of these estimates 
to the exact values as functions of increasing rib, for M S 
{100,1000,10000}. This is shown in Figs. 4-6. 



FIG. 5: The estimated skewness /?3,jv,est, compared with the 
exact skewness /Is.jv of the sampled distribution, as a function 
of rib and M. 



FIG. 6: The estimated excess kurtosis /34,jv,est — 3, compared 
with the exact value — 3 of the sampled distribution, 
as a function of rib and M. For aesthetic purposes, we have 
omitted the large fluctuations of the M = 100 (red) curve for 
rib < 80. 


All of three figures reveal a consistent trend: by in¬ 
creasing the number of bins Ub and/or the number of bin- 
nings M, we can extract with good accuracy the variance, 
the skewness and the excess kurtosis of an event-by-event 
distribution in the observable X. Improving the accuracy 
of this extraction comes with the numerical expense asso¬ 
ciated with increasing rib and/or M, exacerbated by the 
fact that for increasing rib one is decreasing the number 
of events-per-bin n = Ngy/rib, resulting in larger finite- 
number statistical fluctuations of the bin averages. For 
a finite size N of the total ensemble of measured events, 
there is therefore a maximal amount of information about 
the underlying A-distribution that can be extracted from 
the data, even in the limit of infinite computational re¬ 
sources. Experimental HBT analyses typically are based 
on event samples of size 10® — 10^ rather than 10"^. The 
first three moments of the HBT radii distributions should 
thus be accessible with statistical precision of better than 
1%, which is significantly below the typical systematic 
uncertainties associated with HBT measurements that 
are unrelated to the method discussed here. At this level 
of precision, these three moments may already be able to 
provide insights into the nature of the physics lying at 
the origin of the HBT radii fluctuations and help to fur¬ 
ther constrain the material properties of the quark-gluon 
plasma created in relativistic heavy-ion collisions [2]. 


C. Realistic application 

We finally show how the combination of methods 
demonstrated in the previous two subsections can be used 
to obtain physically interesting results. For the same 
5,000 hydrodynamic events considered in Sec. IV, we use 
an azimuthally averaged, ensemble-averaged correlation 
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FIG. 7: For A^ev = ®00 and rib = 2, the left hand panels show 
the estimates for the relative variance for (red), (blue), 
and R 4 (green) for M — 100, while the right hand panels show 
the same resnlts for M = 10,000. Similarly, the upper row 
shows our estimated results for only fib = 2, while the bottom 
row shows the same results for hb = 10. The solid lines are 
the exact results for the ensemble of 5,000 events, while the 
dashed lines represent the average estimate we obtain after 
100 random iterations of our method, and the shaded bands 
indicate the resulting standard deviation of our estimate. 


function method to estimate the coefficient of variation 
(or relative variance) of the event-by-event radii. 

The potential theoretical significance of the relative vari¬ 
ance in heavy-ion collisions has been recently explored in 
[2]. Since the HBT radii are known theoretically on an 
event-by-event basis, we can compute this quantity ex¬ 
actly from Eqs. (33) and (39). We now demonstrate how 
these exact quantities may be reliably estimated by the 
combination of mean estimation and variance estimation 
that we have presented above. 

In detail, the combined algorithm works in the follow¬ 
ing way. We begin with the same ensemble of TVev = 
5,000 events that was already used in Sec. IV, and go 
through M iterations of binning them into bins. For 
each bin, we estimate the arithmetic bin average by 
SMftdividing each bin into ffi sub-bins, and using the 
method presented in Sec. V. Once the arithmetic aver¬ 
age for each of the nt bins is known, we use the method 
presented in Sec. VI to estimate the variance afj. 

We show the results of this combined procedure in 
Fig. 7. Since each of the M binning iterations requires 
a random partition of events into Ub bins, our esti¬ 
mates will tend to exhibit some variability, particularly if 
M is small (say, of order 100). To quantify the resulting 
uncertainty of the estimates shown in Fig. 7, we compute 
the estimates 100 different times, and then compute the 
mean and variance of these estimates. The mean esti¬ 
mates are shown as dashed lines in Fig. 7, and the shaded 
bands represent the ±lcr variability of our estimate. 

For fib = 2, the mean estimate clearly disagrees with 
the exact result for Kt > 1 GeV. This bias is reduced 


by increasing the number of sub-bins: taking fib = 10 
reduces the bias of our estimation procedure almost to 
zero. Similarly, with M = 100, the variability of our esti¬ 
mation procedure is noticeable, but for M = 10,000, the 
widths of the shaded bands are almost negligible. Thus, 
we see that increasing the number of sub-bins fib (top 
row vs. bottom row) results in the ability to decrease 
the bias in our estimate of the exact result, while in¬ 
creasing the total number of binning iterations M (left 
column vs. right column) reduces the overall variability 
of our estimate. For Ub = 2, hf, = 10 and M — 10,000, 
the statistical uncertainty of our estimation procedure ef¬ 
fectively disappears, and the methods presented in this 
paper provide us with a reliable way of accessing statis¬ 
tical moments of the ensemble distribution of the HBT 
radii. 


IX. OUTLOOK 


In this paper, we presented a method for extracting 
from experimental data properties of the event-by-event 
distributions of HBT radii that characterize heavy-ion 
collisions, by estimating their central moments. This 
method allows to overcome the current restriction of pub¬ 
lished HBT results to ensembled-averaged experimental 
data, and opens a way to experimentally access valu¬ 
able information contained in the event-by-event fluctu¬ 
ations of HBT radii that complements analogous infor¬ 
mation from fluctuations of the momentum spectra and 
their associated anisotropic flow coefficients. The pro¬ 
posed method works for both azimuthally averaged and 
azimuthally sensitive HBT analyses. It requires large, 
but not exorbitant event statistics. 

With this work, we promote the area of HBT inter¬ 
ferometry of heavy-ion collisions to a new level that per¬ 
mits the systematic investigation of the statistical proper¬ 
ties of event-by-event fluctuations of interferometric sig¬ 
natures and thereby advances this subfield to a similar 
level of sophistication as established over the last decade 
for momentum-space heavy-ion observables. We expect 
event-by-event HBT analyses to bear similarly rich fruit 
as has been recently harvested from studying event-by- 
event fluctuations of multiplicities and collective flow sig¬ 
natures. 
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Appendix A: Derivation of the variance estimator 

In this appendix, we prove that (42) reduces to (39) when M is taken to be the total number nt) of possible 

distinct ways to sort N events into rib bins of n = N/rib events each. To do this, we first expand the righthand side 
of (39): 


^o,N - ^ 


^ l//^\l/^\ \ 1 -^ N N 

^ fee? fee. S 


i—1 


N{N-l) 


We see that this expression consists of “quadratic terms” (the first sum) and “cross terms” (the second sum). The 
maximal number of independent ways to distribute the total ensemble of N events into rib bins of size n = N/rib is 


M^^^{N,nb) = 


N\ fN- 


n / \ n 


2n\ _ Nl _ m 
n) (n!)”'> ((A/n;,)!)"'’ 


For M = Mmax, we may write the righthand side of (42) as 




nb{nb-l)M„ 




j = l k=l 


A^max(n& —1) 


f" ((0)5 - (o)t). 


where 

•^^max -^^max -^^max -^^max 

E (‘5>h = E <«’>?.. = ■ ■ ■ = E ^ E («’>? ■ (Ai) 

i=i i=i i=i i=i 

This equality holds by definition when M^ax exhausts all binning possibilities, because the different sums in (A4) 
only differ by permutations of their summands. 

Introducing the notation to represent the fcth event observable in the jth iteration,® we can expand both of 
the terms in (A3) as follows, to bring them into the form (Al). Using (A4), the first term in (A3) becomes 


TVT AT ^'^m.ax / -< ^ 

_^_V io?- = _ - _V ( - 

A^max(?^6—1) ■“ ^ Af^lax(?^6—1) ■“ \ ^ ffe ^ 


^max 


_E_V V V 

n®M_(n,-I) - 


Mmax r n 


n^Mniax(?T'fc 1) 


lyE EK'T + 2EE^i''^^i'^ 


Afiyiax(^b 1) 


Ol ^ 02+202 


oi is a degeneracy factor which counts the number of times any given event falls into the first bin; with our choice of 
Afmax) this is just Mffjax/u;,. Similarly, 02 is a degeneracy factor which counts the number of iterations for which two 
different event both fall into the first bin; it is given by 

fN-2\fN-n\ (2n\ Mi„axn(n-I) 

“^^(„-2)( „ )--q.,)= jv(jv-i) <“> 


^ Of course, since each iteration can be thought of as a partition 
into Til) bins of a random permutation of all N events in the 
ensemble, the kth event observable in the jth iteration will gen¬ 
erally not be the same as kth event observable of the j + 1st 
iteration. 
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The second term in (A3) is less subtle: 


, N N 


A^max(?T-6 “ 1) 
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Combining these results, we obtain 
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(AS) 


which is identical to (Al). This establishes that cro^jv^est fro m (42) reduces to cfq j^ from (39) when M = 


Appendix B: Two derivations of 


particle and two-particle spectra is written 


In this appendix, we present two different ways of un¬ 
derstanding the relationship between the event-by-event 
HBT radii and the experimentally measured, ensemble- 
averaged radii, . First, we show how the relationship 
(10) arises from the combination of (1) and (8). We then 
show that precisely the same result holds for the radii ex¬ 
tracted from (14) via the source-variances method, mo¬ 
tivating us to treat (14) as the most physically accurate 
way of accounting for the contribution of event-by-event 
fluctuations to the experimental correlator (10). 

We begin by introducing the convenient shorthand 


(C')ev = 


(A^1.2) 

{Ni) (iVa) 
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d^N^k) 
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£ = 1 , 2 , 


(B2) 


where k labels a particular event in the ensemble. In this 
notation, the correlation function (1) for the fcth event is 
written 


C(k) = 


jy(k) 


(fc) Ar(fe) 


ivr'w 


= 1-1- exp 



Eij Qi Qj 


(B3) 


where we defined 


Wk = 


N, 

ev 1 


(k) j\^ik) 


(e 


iVev 

k^l 



N^k)Nik) 


(B5) 


The last step follows in the smoothness approximation 
[17] if we additionally define N = (Ni) = {N 2 ). 

Although the single event correlators (B3) reflect the 
traditionally imposed normalizations 


lim {q, k) = 2, Jim (g, k) = l, (B6) 

g—>-0 q—¥oo 


the ensemble-averaged correlator defined by Eq. (B4) 
leads to a different overall normalization: 


where we have suppressed the functional momentum de¬ 
pendence for the sake of clarity. On the other hand, 
the correlator defined from the ensemble-averaged one- 


(Nk ^ (Nk 

jTo (.-.^)=2V’ ir. (.1 ^) = V- 

(B7) 
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Consequently, we fit (B4) to a Gaussian form with an 
additional normalization factor: 


(C')g^ = Co 1 + exp - 






(B8) 


where we have to set Cq = (-^^) in order to satisfy 
(B7). Then, equating (B8) with (B4) and requiring their 
curvatures (with respect to q) to be the same as (f —)■ 0 
leads immediately to the expression 


and 


{C)^,{q,K) = 


1 + 


1 + 

1 + 


1/ d‘ixS{x,K)\^ 

d^xe^'^-So{x,Kf 
\J d‘^xSoix,K)\^ 

(C-1) (B14) 


_/V2 (*J> ]\[2 > 


(B9) 


which is simply Eq. (10). 

We can derive this same result by a slightly different 
treatment. Instead of directly requiring the ensemble- 
averaged correlator to satisfy certain normalization con¬ 
straints, we now show how to obtain the source-variance 
HBT radii from the physical correlation function (14) 
which correctly incorporates the complete effects of 
event-by-event fluctuations. By utilizing the source- 
variance method, we choose to approximate the corre¬ 
lation function as a Gaussian in q, which means the cor¬ 
responding radii may be extracted as above by simply 
computing the curvature of the correlator as g—)■ 0: 




(C), 


2((C')ev-l) dq^dqj 


(BIO) 


q-^O 


Here, the factor of (C)^^ — 1 in the denominator of this 
expression account for the eventwise multiplicity fluctua¬ 
tions of A^. To see this, consider an ensemble of emission 
functions which differ only in overall normalization N: 


Sn{x,K)=NSo{x,K), (Bll) 

where So{x,K) has the same Gaussian form for every 
event in the ensemble. If the fluctuations of normaliza¬ 
tion are governed by a probability distribution P{N), 
then we can write the ensemble-averaged emission func¬ 
tion S{x, K) as 

S{x,K) = {Sn{x,K))^^ 

= So{x,K)(^J dNP{N)N^ 

= NSoix,K), (B12) 

where N is the average normalization for the ensemble. 
Then, using Eqs. (12) and (14), it is straightforward to 
show that 


C{q,K) = l + 


fd^xe^‘^-^So(x,K) 

fd^xSo(x,K)j^ 


(B13) 


This result demonstrates that, although (C}^^ and C 
both have the same dependence on q and therefore the 
same HBT radii, Pfj = , the curvatures of the two 

correlators do not agree as (f —)■ 0: rather, they differ 
by an additional factor of which is equal to 

(G)g^ ((f —>■ 0) — 1. Gonsequently, we define the radii in 
terms of the appropriately normalized curvature of the 
correlator at the origin.Using Eq. (BIO), we can then 
write 


AT,, 


^%3) 




/ d^xd^x' 


^ K)Skix', K) 


q-^0 


—^ / d^xd^x'Skix,K)Skix',K) 

^ /ev J 

X [ (Xi-Pit) - {x'^-Pit') ] [ (Xj-Pjt) - {x'j-Pjt') ] 

1 j\i2 , 

7^ ^ (((a;i-/3*t) {xj-Pjt))g^ 

— {{Xi — Pit)) {i^j~ Pj^)) Sk} 

(B15) 




{N^)e 


where, in the last step, we have simplified our result by 
means of Eq. (6). This result is, again, seen to be equiv¬ 
alent to (10). 

This means that a consistent extraction of the 
ensemble-averaged HBT radii should account for the fact 
that more pairs of identical particles will come from 
events with larger total charged multiplicities dN^^/dq, 
leading to the HBT radii for these events being more 
heavily represented in the final, ensemble-averaged radii. 
Consequently, the true average radii are in fact a weighted 
average over the event-by-event HBT radii, where the 
weighting factor of / (iV^) represents the fraction of 
all pion pairs used in the construction of (8) that come 
from an event with overall normalization N. 
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